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Abstract. Motivated by the role that spectral properties play for the dynamical 
evolution of a quantum many-body system, we investigate the level spacing statistics 
of the extended Bose-Hubbard model. In particular, we focus on the distribution of 
the ratio of adjacent level spacings, useful at large interaction, to distinguish between 
chaotic and non-chaotic regimes. After revisiting the bare Bose-Hubbard model, we 
study the effect of two different perturbations: next-nearest neighbor hopping and 
nearest-neighbor interaction. The system size dependence is investigated together 
with the effect of the proximity to integrable points or lines. Lastly, we discuss the 
consequences of a cutoff in the number of onsite bosons onto the level statistics. 
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1. Introduction 

In classical systems, the chaotic nature of a Hamiltonian plays a crucial role in 
determining the dynamics. Chaotic systems explore a large area of their phase space, 
whereas non chaotic ones can be trapped in certain atypical subspaces p. Thus, chaos 
plays a considerable role in explaining thermalization and ergodicity. For quantum 
systems the situation is less clear. In principle, the quantum dynamics is linear, since 
it follows Schrodinger equation, and therefore the notion of chaos is not so well defined. 
However, in quantum systems for which a classical counterpart exists, one finds that the 
dynamics is very different depending on whether the corresponding classical motion is 
chaotic or regular [1] . For generic quantum systems, chaos is also believed to be essential 
for thermalization [21 [3] and delocalization in Fock space [1] . The recent realization of 
closed quantum systems out of equilibrium by strongly correlated cold atoms clouds 
[5] have lead to a renewal of interest in the chaotic properties of many-body quantum 
systems and their relation to thermalization [6]. 

For chaotic quantum systems it has been conjectured that their spectra show 
universal features which are related to the theory of random matrices (Tj [HI [HI [TO] . 
To quantify the spectral properties of a given Hamiltonian, a natural quantity to look 
at is the gap between adjacent many-body levels Sn = E^+i — En, where {En} is the 
list of eigenvalues in ascending order. The general symmetries of the Hamiltonian, like 
time-reversal and half-integer spin rotational invariance, determine the random matrix 
ensemble (among orthogonal, unitary and symplectic ensemble) to which it belongs [?]. 
A spin-less time-reversal invariant Hamiltonian such as the Bose-Hubbard model (for 
generic parameters) should have similar universal features as the Gaussian orthogonal 
ensemble (GOE). For example the adjacent level-spacing distribution Pa{S) is predicted 
to take the Wigner-Dyson form 



where A is the mean level spacing. On the contrary for so-called integrable models, in 
which the properties of the system are determined by an extensive number of conserved 
quantities, the level-spacings should exhibit the following Poissonian distribution 



The relation between the spectral properties of a quantum system and the random 
matrix ensembles has been demonstrated numerically even for strongly correlated many 
body systems without classical counterpart. In particular a GOE like behavior was 
pointed out for the non- integrable two dimensional t-J model [11]. In one dimensional 
models, it was checked that at and close to the integrable points, the statistics was 
Poissonian-like [12] . In one dimension, non-diffractive models are integrable on rigorous 
grounds and solvable using the Bethe-ansatz which provides all eigenstates [13] . Several 
other systems have been discussed since then [HI [IDl [iSl [1]. Therefore, analyzing the 
spectral properties of a system provides a phenomenological approach to investigate the 
chaotic nature of the quantum Hamiltonian. 




(1) 




(2) 
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In this work we consider the properties of the spectrum of the Bose- Hubbard model. 
This model is a paradigmatic strongly correlated many-body system where interactions 
of amplitude U compete with the kinetic energy favored by the hopping J. It appeared 
in several contexts of condensed matter theory and regained a lot of interest since it 
has been reahzed in quantum gases confined to artificial lattice structures [5]. Its out- 
of-equilibrium dynamics and in particular the question of thermalization following a 
quantum quench have also been studied numerically [HI |T71 [HI [191 |2Q] . As perturbations 
to the Bose-Hubbard Hamiltonian are experimentally relevant in different regimes [2T] . 
it is essential to address the sensitivity of the spectral features to extra terms. We 
consider the effect of two different perturbations and their consequences on the level 
statistics. We focus on the model with a next-nearest-neighbor hopping with amplitude 
J2 as well as a nearest-neighbor interaction with amplitude V . The Hamiltonian under 
study then reads: 

L L 

i=i j=i 
U ^ L 

where 6j and bj are the bosonic creation and annihilation operators, and hj = h^-b^ 
the number operators on site j and L the number of sites in the chain. Periodic 
boundary conditions are used to take advantage of translational invariance. We analyze 
the spectrum in the subspace spanned by reflection symmetric states with zero total 
momentum using full exact diagonalization. No cutoff in the onsite occupation M is 
assumed, i.e. M = N, and unit filling N/L = 1 is taken if not stated otherwise (A^ 
is the total number of bosons). The symmetrized Hilbert space has dimensions up to 
56, 822 for N = L = 12. The considered extended Bose-Hubbard model is integrable (in 
a broad sense) only in the two limiting cases of ?7 = = (free bosons) and J = J2 = 
(atomic limit). 

In a first analysis on small lattices, Kolovsky and Buchleitner [22] discussed the 
chaotic behavior of the bare Bose-Hubbard model (V = J2 = 0). They investigated the 
level statistics and the Shannon entropy of the eigenstates with respect to the basis of 
the two integrable limits as a probe of the eigenstates delocalization (another chaotic 
feature). A pronounced similarity of the level statistics with GOE and a maximal 
Shannon entropy has been found when the hopping amplitude and the interaction 
strength are of the same order. Furthermore, it has been noticed that the low and 
very high energy parts of the spectrum seem to display features of integrable spectra. A 
chaotic trimeric version of the Bose-Hubbard model has also been analyzed [231 121] and 
a comparison to the semi-classical approximation has been carried out when N/L S> 1. 
Another case that has been studied corresponds to restricting the maximum onsite 
number of bosons M to one. In this case, the Bose-Hubbard model boils down to a 
hard-core boson model (XXZ) which has the particularity of being integrable provided 
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Figure 1. The level spacings distribution ioi U/J — 2 and U/J — 40 for a system 
of length L = 12 are shown (symbols). The dash-dotted, green curve is the Poisson 
distribution and the dashed, red curve the Wigner Dyson distribution of a purely GOE 
like ensemble. The different symbols correspond to different unfolding procedures. For 
the curves labeled by the number of energy bins, the unfolding is performed locally 
over energy regions, which divides the total energy range covered by the eigenvalues 
into equally spaced energy bins. The smooth part of the level staircase is fitted by 
a polynomial. Additionally an unfolding within the energy bands is performed for 
U/J = 40 using the function A2 + {Al - A2)/{1 + exp{{E - Eo)/AE)) where Al, A2, 
Ef) and AE are fitting parameters. 



J2 = 0. The level statistics of this hard-core model, perturbed by a non-integrable 
operator, has been studied recently |25] . While the two models are physically equivalent 
in the low-energy physics provided that U is large enough, the high-energy spectra are 
very different even in the large-f/ limit. 

In this paper, we analyze the level statistics of the extended Bose-Hubbard model 
in detail, using system sizes up to L = = 12 sites. All data are restricted to the zero- 
momentum and reflection symmetric sector (the one of the ground-state). In section 
|2] we study both the level spacings distribution and the distribution of the ratio of 
consecutive level spacings. The latter is a very useful observable, in particular close 
to the atomic limit. This allows us to address the energy and finite size dependences. 
Furthermore, we investigate in section |3] the influence of different perturbations on the 
spectrum of the Bose-Hubbard model. We consider a next-nearest neighbor hopping 
and a nearest neighbor interaction and focus is put on how the chaotic features evolve 
away from integrable points (lines) in the parameter space. Lastly, we discuss the effect 
of the onsite occupation cutoff M on the level statistics. 

2. Properties of the bare Bose-Hubbard model 

We first focus on the bare Bose-Hubbard model (J2 = V = 0) [22j. In Fig. [H typical level 
spacings distributions of the unfolded spectrum are shown for two different interaction 
strengths. In order to remove the dependence on the system specific mean level density, 
the 'unfolding' procedure consists in renormalizing the level spacings by using a suitable 
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Figure 2. The many-body density of states p in the zero- momentum and reflection 
symmetric sector and for different interaction strengths. Eq is the ground-state energy. 
The evolution from a smooth spectrum at low interaction strength to the appearance of 
well separated energy bands at large interaction strength is shown. The smooth curves 
are obtained broadening the (5-peaks by Gaussians {p{E) — J2n9^i-^ " ^n) with gA 
a Gaussian of width A). However, note that even for large interaction strengths the 
width of the Gaussians is much smaller than the width of the energy bands. 

fit for the smooth part (see Ref. |10j for details of the procedure). The left panel of Fig. [1] 
shows that the distribution for U / J = 2 closely follows the Wigner-Dyson distribution 
of the GOE. In particular, the distribution vanishes for small level spacings, a typical 
signature of level repulsion associated with avoided level crossings. When approaching 
the integrable points the distribution deviates from Wigner-Dyson. This is shown for 
strong interaction f// J = 40 in the right panel of Fig. [T|, but occurs also when lowering 
the interaction. The tail of the distribution is close to the exponential tail of the Poisson 
distribution; for small level spacings there is a significant enhancement compared to the 
GOE distribution but some level repulsion persists. We observe that increasing the 
size of the system (not shown, see also the discussion here after) tends to increase the 
similarity to the Wigner-Dyson distribution. In particular, the repulsion at small level 
spacings becomes more and more pronounced. However, for currently accessible system 
sizes, it is not clear whether the distribution very close to the integrable points will 
converge to the Wigner-Dyson one. 

The study of systems with large interaction strength is involved due to the 
appearance of a band structure in the density of states (see Fig. [2]). The spectrum 
evolves from a smooth broad spectrum for small interaction strength to a series of 
narrow energy bands separated by U for large interaction strengtli||. This would make 
the unfolding complicated because it is difficult to separate the spectrum into a smooth 
and a fiuctuating part. To obtain reliable results we performed different unfolding 
procedures fitting locally different parts of the spectrum. If the ranges over which the 



X However, note that the mean level spacing is still much smaller than the width of the energy bands. 



Statistical properties of the spectrum of the extended Bose-Huhbard model 



6 




Figure 3. Left panel: semicircular density of state for an GOE ensemble of random 
matrices of size 5, 000^ averaged over 10, 000 realizations. In addition, the value of the 
ratio of consecutive level spacings, averaged over the different realizations, is given as 
a function of energy, showing that it does not depend on the density of states. The 
shaded regions are the regions over which the distributions shown in the central and 
right panel are taken (they both contain 200 states). Central panel: Distribution of 
the ratio of consecutive level spacings taken over a region at the edge, in the bulk and 
of the whole spectrum: the distribution does not depend on the variation of the local 
density of states. Right panel: Distribution of the level spacing in the same regions 
without unfolding the spectrum showing the strong dependence on p{E). 



fits are performed are chosen in a suitable way, the same general form of the spectra is 
recovered, even though small discrepancies can occur. Note that the complicated form 
of the density of states also makes a study of longer range correlations of the spectrum 
involved. 

To avoid this complication, we continue our study using another measure which 
has the advantage of not depending on the unfolding procedure: the ratio of 
consecutive gaps between adjacent levels [26]. This quantity is defined by r„ = 
min(5„, max((5„, As long as the density of states does not vary on the 

scale of the mean level spacing, the trivial dependence on the smooth part drops out 
and there is no need for unfolding. This is exemplified in Fig. |3l where we compare 
for a GOE random matrix ensemble the distribution of the level spacing and of the 
ratio of the consecutive level spacings taken over two different energy ranges without 
performing the unfolding procedure. The first energy range lies at the boundary of the 
spectrum where the density of states varies considerably. In contrast the second energy 
range is situated in the center of the spectrum, where the density of states exhibits 
only slow changes. The GOE distribution is computed numerically as in Ref. [26] using 
the averaged results of 10, 000 samples of random matrices of size 5, 000^. It is clearly 
seen that the ratio of consecutive level spacings does not depend on the region used, 
whereas the level spacing statistics is different for the two chosen regions due to its 
strong dependence on the smooth part of the spectrum. The local ratio of consecutive 
level spacings averaged over the GOE ensemble is also independent of the density of 




Figure 4. Typical distribution of the ratio of adjacent level spacings r. Left panel: 
Different interaction strengths for L = 12 are compared with Poisson and GOE 
distributions. The proximity to the chaotic regime is seen when [/ ~ J. While for 
U ^ J and U ^ J the curves approach the Poissonian tail. Right panel: system size 
dependence for U / J = 40. 



states (left panel of Fig. [3]). On the basis of these results, which show that r is a useful 
random variable to analyze the spectral statistics, we continue our study focusing on 
the ratio of consecutive gaps. 

We give in the left panel of Fig. H] typical distributions P(r) obtained for the Bose- 
Hubbard model at different interaction strength. As expected, the P(r) distribution 
depends on the chaoticity of the Hamiltonian: the distribution for a Poissonian spectrum 
reads P(r) = 2/(1 + r)^ while it can been computed numerically for the GOE ensemble 
(Fig. [3] and Fig. H]). As for the level spacing distribution, the maximum resemblance 
with GOE is observed when f/ ~ J. In the proximity of the two integrable limits f/ = 
and J = the distributions approach the Poisson prediction, particularly on the tail. 
We checked for a wide range of values of U that the distribution of these ratios is in good 
agreement with the level spacing distribution using different local unfolding procedures. 

In order to have a systematic tool to probe the proximity either to Poisson or 
GOE, we use the results on the mean value (r) which is (r)p = 2 ln2 — 1 ^ 0.386 for 
Poisson and {r)GOE ~ 0.53(1) for GOE. We thus expect that this averaged ratio should 
display a maximum for an interaction strength around U ^ J. In Fig. |5] this ratio is 
given for different interaction strengths and system lengths. At intermediate interaction 
strength (0.3 < fZ/J < 4) we see that (r) gets very close to the GOE prediction. Even 
though there is a small system size dependence left, the values for the longer system 
sizes considered are quite close to the expected value. For small and large values of the 
interaction strength we see that the behavior is different. For these regimes the values of 
(r) lie in between the ones expected for the GOE and the Poisson ensemble and for most 



Statistical properties of the spectrum of the extended Bose-Huhbard model 



8 



0.55 



0.45 



0.35 




Figure 5. Evolution of the average ratio of consecutive level spacings (r) with 
the interaction strength U/J. For the Poisson distribution the average value of 
(r) p = 0.386 and for the GOE ensemble {r)GOE = 0.53. The average is taken for 
the full spectrum, without cutting the low and the high energy part. A clear tendency 
of larger system sizes versus the GOE ensemble value is seen. 
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Figure 6. Scaling of the ratio of consecutive level spacings (r) with the inverse system 
size for different interaction strength. 



of these values a strong system size dependence is still evident. Typically, the trend for 
longer system sizes goes towards {r)GOE, as shown in Fig. [6]for some chosen values of the 
interaction. This suggests that away from the integrable points some critical length scale 
(or particle number) exists above which the system shows a level statistics which is very 
close to the one of the GOE ensemble. Our results further suggest that this length scale 
possibly grows in the proximity of the integrable points. For instance, when U/J = 40, 
the distributions hardly evolve with the system size (see Fig. H] right panel). We expect 
that for large enough sizes the properties of the spectrum might be well described by a 
GOE. However, the question whether or not a finite deviation of the parameters from the 
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Figure 7. Average ratio of consecutive level spacings taken in different energy 
intervals of equal size for different interaction strength U/J — 0.1,2,40 and different 
system length L = 11, 12. The number of bins taken for the average is 50. 



integrable limit is necessary to obtain GOE like characteristics cannot be conclusively 
answerecl§]. Still, the obtained results can be compared with other scenarios on the effect 
of the interaction in the Bose-Hubbard model. For instance, Cassidy et al suggested [27] 
that there could be an interaction threshold in the Bose-Hubbard model for the chaotic 
behavior to develop, based on calculations valid in the semi-classical limit N/L ^ 1 
supplemented by a mean-field calculation. Extrapolating their results to the N = L 
limit, the threshold would be f//J ^ 0.5. In contrast, our results demonstrate that for 
low filling n = 1 even at U/Jc:^ 0.1, the level statistics features have a strong tendency 
towards a chaotic behavior and no signature of a threshold is found. 

Up to now we have considered the properties of the whole spectrum. However 
the question arises how these properties change within different energy ranges of the 
spectrum. To illustrate this, we show in Fig. [7] the average ratio taken over different 
ranges of the spectrunIJi]. We must notice that contrary to the GOE benchmark of Fig. [3l 
the Hamiltonian is deterministic and no sampling can smoothen the curves: we were 
able to get good statistics only by reaching large enough systems sizes. For small values 
of the interaction strength (f// J = 0.1 and f// J = 2 in Fig. [7]) the average ratio slightly 
depends on the energy and one can observe that the bulk of the spectrum displays 
the GOE prediction while edges show some deviations with strong fluctuations. These 
strong fluctuations may be attributed either to the fact that the small density of states 
induces bad statistics or to the fact that the physics at the edges (in particular close 
to the ground-state) display different statistics than for the high energy excited states 

§ Notice that the scaling analysis of the level statistics have some strong numerical limitations. Indeed, 
as the width of the spectrum scales as N'^ while the number of states scales exponentially with N, we 
may expect the minimal level spacing to reach the numerical accuracy of full diagonalization at some 
relatively small system size L. However, these system sizes are longer than the here considered system 
sizes. 

II We checked that the ratio exhibits the same features as the level spacing statistics. 
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Figure 8. The Shannon entropy in the real-space symmetrized basis of exact 
diagonalization is shown for different interaction strength. An evolution from a smooth 
behavior for small interaction strength towards a very fluctuating behavior for large 
interaction strength can be seen. 
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Figure 9. The finite system dependence of the average value of the Shannon entropy 
is shown. A clear trend towards larger values is found. 



in the bulk. A clear maximum exists in the central region. Increasing the system size 
for U/J = 0.1 makes the bulk value closer to the GOE prediction. For intermediate 
values of the interaction strength U / J = 2 (central panel Fig. [7]) the properties of the 
spectrum do not change much in different energy regions and are close to GOE. For 
large values of the interaction strength (lower panel in Fig. [7]) a band structure develops 
in the energy spectrum. The ratio shows stronger fluctuation^ and a slight drop from 
a more Wigner-Dyson like value towards a more Poisson like value. 

In addition to the distributions discussed above, the Shannon entropy can give 
interesting information about the chaoticity of the system |22]. The Shannon entropy 
is defined by Sn = (— X^m kmH/l'^dimif, where are the coefficients of the 



^ The ratio and its fluctuations also depend more strongly on the energy interval used. However we 
checked that the main trend remains for typical energy ranges taken. 
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nth eigenstate decomposed onto the mth basis state of a chosen basis. The entropy 
measures the delocahzation of a wave-vector with respect to a chosen basis: with our 
definition, it reaches 1 for a fully delocalized state. Here we measure the entropy with 
respect to the symmetrized real-space basis of exact diagonalization. In Fig. [8] the 
behavior of the entropy is shown for different interaction strengths. A clear evolution 
from a smooth to a sharp distribution around high values at small interaction strength 
towards a strongly fluctuating distribution at low values is evident. This shows that 
for small interaction strength almost all eigenvectors are delocalized as expected in the 
non-interacting regime where the Hamiltonian is diagonal in the momentum space. At 
intermediate interaction strength, the edges display less delocalized features than in the 
bulk. The most interesting behavior is exhibited at large interaction strength where the 
degree of localization fluctuates strongly between different eigenvectors within a Mott 
lobe. Let us point out that this finding is similar to the behavior of some observables 
and weights calculated in these eigenstates and that was identified to be the reason 
of non-thermalization after a quench on such finite size systems [iHl [iHl ED]. Finally, 
let us comment on the evolution of the above results increasing the system size. We 
find that even though Fig. [9] displays a trend towards delocahzation, the behavior in 
the thermodynamic limit is particularly hard to access close to the infinite-?/ integrable 
point. 

3. Perturbing the Bose-Hubbard model 

We now turn to the effect of the J2 and V perturbations (separately) that are expected 
to help breaking integrability at the J = and [/ = integrable points, respectively. As 
there are three parameters ranging from zero to infinity, we will fold the parameter space 
using two representations (see Fig. [10] for an example). The first one is to introduce 
the function f{x,y) = x/{x + y) and to use the following definition: u = f{U, J) and 
j2 = f{J2, J) when J 7^ 0; M = /(f/, J2) and j2 = /(^2, U) ii J = 0, with the additional 
point u = i2 = 1 when U = J2- Such a folding is useful to restrict the considered 
parameters onto a finite interval, and enables one to easily deduce the parameters for 
a given point. One disadvantage of this folding are discontinuities arising from the 
infinities on the j2 = 1 and u = 1 lines. A continuous way to draw the data is to use a 
ternary plolQ. However, it is more difficult to find in such a plot the original parameters. 
Therefore, we use both ways to present our results. 

Influence of a next-nearest neighbor hopping - If the next-nearest neighbor hopping 
J2 is switched on, the behavior of the spectrum changes. A summary of the effect of 
J2 for L = 11 is presented in Fig. [TOl For small and intermediate interaction strength, 
the additional finite value of J2 drives the small systems closer to the Poisson behavior. 

+ Formally speaking, a point of the diagram corresponds to percentage of each parameter, i.e. a 
triplet (% J, % J2, %J7). In cartesian coordinates with the triplet (100,0,0) at the origin, one has 
x = {U + Jil2)IN and y ^ ^/2,J2/2J\f with M ^ J + J2 + U. 
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Figure 10. Density plots of the evolution with the parameters j2 and u (see 
text) of the averaged ratio of consecutive level spacings (r) {L = 11). We plot 
( (r) — (r) p ) / ( (r) p — {r)GOE), so that blue (resp. red) corresponds to Poissonian (GOE) 
distribution. Right panel: ternary plot of the same data showing the integrable lines 
J — J2 and the isolated integrable point {U corner) . 

For large interaction strength, J2 helps to drive the system away of the integrable point 
J = 0, before at very large 72 it again turns Poisson like due to the attraction of the 
U = J = integrable point. The values of J2 can thereby be much smaller than 
the actual interaction strength and still have a considerable influence. In the ternary 
diagram the almost symmetric behavior of the system with respect to the diagonal 
J = J2 is nicely visible. This means that the next-nearest neighbor hopping J2 has a 
similar effect than the nearest neighbor hopping J. In order to discuss the influence of 
longer system sizes, we show in Fig. [11] the dependence of the average ratio (r) on j2 
for different system sizes at chosen values of the interaction strength. If the behavior is 
already close to the GOE one, finite size effects are typically very small (Fig. ITT] central 
panel). In contrast, if the system is not in the GOE regime, the finite size effects become 
more pronounced. However, the larger system sizes show a clear tendency towards the 
GOE behavior. Lastly, we have checked that the same qualitative features are displayed 
in the Shannon entropy when J2 is taken into account. For instance. Fig. [12] shows that 
starting from a large U/J and increasing J2 tends to make the wave-functions delocalize 
in the symmetrized real-space basis of exact diagonalization. In particular, for the point 
with J = J2, the perturbation clearly favors delocalization. When J2/U increases, the 
delocalization with respect to the real space basis becomes more and more pronounced 
due to the dominating kinetic term. As for the case of the pure Bose-Hubbard model, 
a strong dependence on the energy is observed. A maximum in the Shannon entropy 
is found for intermediate energies, whereas for low and in particular high energy the 
Shannon entropy drops quickly and exhibits typically more fluctuations. 
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Figure 11. Evolution of the ratio of consecutive level spacings (r) with nn hopping 
j2 for different interaction strengths. For the Poisson distribution the average value of 
(r)p = 0.386 and for the GOE ensemble {r)GOE — 0.53 are displayed. The average is 
taken for the full spectrum, without cutting the low and the high energy part (circles) . 
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Figure 12. Effect of the J2 perturbation on the Shannon entropy in the symmetrized 
real-space basis. 



Influence of the nearest-neighbor interaction - In the following paragraph we discuss 
the influence of nearest-neighbor interactions which is summarized in Fig. [131 The same 
representations of the parameters space are taken, replacing J2 with V and j2 with v. 
The behavior is qualitatively very different from the J2 perturbation: the integrable 
line in the ternary plot is now the U — V line. In both representations, the data are 
nearly symmetrical with respect to the diagonal U = V . Thus, the two interacting terms 
act in a similar way in terms of level statistics. For small onsite interaction strength 
f/, a finite nearest neighbor interaction V enhances the trend of the ratio towards its 
GOE value. However-surprisingly at first sight-at large onsite interaction strength [/, 
a small finite value of V induces a trend towards the GOE like behavior and only if 
V is larger than the onsite interaction U the value of the ratio drops drastically to the 
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Figure 13. Same as Fig. [TOl changing J2 with nearest-neighbor interaction V. 



Poisson value. In particular for interactions of the same order of magnitude U ^ V, J 
rapidly drives the system towards GOE. This effect can be made plausible in a simplified 
picture considering the eigenstates of the Hamiltonian at J = 0, which are Fock states. 
However, their order with respect to the energies is very different for both interactions. 
To make this more exphcit consider the state with one particle per site |1) = |1, 1, . . . , 1) 
and the state with two particles every second site |2) = |2, 0, 2, 0, . . . , 2, 0). For a strong 
onsite interaction the state |1) is very low in energy whereas the state |2) lies in the 
upper part of the spectrum. In contrast, for a strong nn-interaction state |2) lies in the 
lower part of the spectrum whereas state |1) lies in the upper part. If both interactions 
are of the same order of magnitude both states lie very close in energy such that the 
small hopping has a large effect on the states. These two energy states are examples of 
the behavior of many of the energy states which become almost degenerate in the limit 
of equal onsite and nn-interaction strength. Therefore, the effect of the hopping as a 
perturbation is expected to be more effective when U ^ V and should help make the 
level statistics GOE-like. 

In Fig. [131 finite size effects are considered. As for the other discussed cases the 
finite size effects are very small if the value of the ratio is already close to GOE. For the 
remaining values, we typically see a trend of the ratio for larger system sizes towards 



Occupation cutoff dependence - In Fig. [15], we show the effect of introducing a cutoff 
in the number of bosons M on each site. For large values of the interaction strength the 
use of a cutoff M > 4 does not seem to change the spectral properties much, both at 
unit and half fillings. A smaller cutoff strongly affects the ratio, pushing it close to the 
GOE prediction. This is due to the suppression of the remaining particle fluctuations 
mixing the true eigenstates which are no longer represented. For intermediate interaction 
where the ratio is close to the GOE value, the effect of the cutoff on the mean ratio is 



GOE. 
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Figure 14. Same as Fig. [Til changing J2 with nearest-neighbor interaction V. 



relatively small. In contrast for small interaction strength the influence of the cutoff is 
very pronounced. Here the introduction of a small cutoff does drive the system away 
from integrability. Even for [7 = 0, using M < N makes the levels statistics close to 
GOE. This effect can be understood by recalling the properties of the eigenstates in 
the limit of weak interaction. These are the momentum eigenstates which comprise 
strong particle fluctuations. If one introduces a cutoff for the number of bosons per 
site these states cannot be represented anymore and start to mix. In other words: the 
local constraint is equivalent to using a projector on the kinetic part which correlates 
the bosons or, equivalently, acts as a complicated effective interaction that turns out 
to display a GOE behavior. A qualitatively similar effect is found in the comparison 
between the ID t-J and Hubbard model. The Hubbard model is integrable while the t-J 
model which is related to the Hubbard model by Gutzwiller's projection is generically 
not [12]. 

In addition, we can discuss earlier results in the literature addressing the question 
of the integrability of the ID Bose-Hubbard model and the effect of multi-occupancies. 
Seminal studies by Choy and Haldane [28| 129] seemed to argue that Bethe-ansatz 
equations yield solutions of Bose-Hubbard-like models but the analysis turned out 
to be invalid |30l |3T]. The authors emphasized that M > 3 was required to give 
rise to non-integrability. Later, Krauth [32] used the Bethe-ansatz wave-function as 
a variational approach for the ground state properties. He found that for f//J < 2, 
the comparison with unbiased quantum Monte-Carlo results was indeed very good. 
This finding supported the fact that the integrable nature of the free bosons gas was 
preserved up to interactions close to the transition point to the Mott insulating phase 
at least for ground state properties. The results of the present study, in which we 
found considering the entire spectrum that the chaotic properties emerge much below 
f//J ~ 2 (u ~ 0.66), stress the difference between the low-energy part of the spectrum 
and high-energy regions: the ground-state and first excitations might have integrable- 
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Figure 15. Evolution of ((r) — {r)p)/{{r)p — {r)GOE) witli tlie maximum number 
of onsite bosons M . The average is taken for the full spectrum, without cutting the 
low and the high energy part (circles). Left panel: incommensurate density where 
hard-core bosons corresponds to M = 1. Right panel: the same at a commensurate 
density N — L. 



like behavior (if the density of quasi-particles is small, they may interact less) while one 
cannot consider a high energy excitation as simply being made of a superposition of 
elementary excitations [11] (a picture which survives high in energy in the Bethe-ansatz 
and in free particles systems). 



4. Conclusion 



To conclude, we presented a study of the characteristic properties of the spectra of the 
extended Bose-Hubbard model. In an intermediate regime of the interaction strength the 
system is in the GOE regime. In contrast for very weak and strong interaction strength 
the analysis suggests an approach toward GOE when increasing system sizes. In most 
parameter regimes this trend towards the GOE was most pronounced in the central 
region of the spectrum. An additional next-neighbor hopping amplitude J2 changes the 
properties of the energy levels in these small systems. It acts similar to the hopping 
amplitude J. For weak interaction, J2 drives the system closer to the Poisson like 
behavior, whereas for large interaction strength it reinforces the GOE like behavior. An 
additional nearest neighbor interaction V has a similar effect on the statistical properties 
of the spectrum as the onsite interaction even though the corresponding eigenstates are 
very differently distributed in energy. Close to the point where the interaction U and 
V become of similar strength even a very small value of J is enough to induce a GOE 
like statistics. Finally we discussed the influence of the introduction of a cutoff for the 
number of bosons per site usually used to render the system numerically tractable. Here 
we see that the cutoff can change the statistics of the spectrum from Poisson like to 
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GOE like, in particular at small interaction strength. 

We see that for all the different regimes considered the changes with increasing 
system size can be divided into two main classes. If the properties of the system are 
already GOE like, increasing the system size only induces small changes. In contrast if 
the value of ratio of consecutive level spacings lies in between the Poisson and the GOE 
value indicating a mixed statistics, finite size effect are considerable. In this regime larger 
system sizes typically drive the system towards the GOE value indicating a GOE like 
behavior in the thermodynamic limit. However, larger sizes would be needed to obtain 
a conclusive result on the question whether there is always a large enough system size to 
reach a GOE behavior for all parameters except the ones corresponding to the integrable 
points (or lines) or if a threshold for the perturbation from the integrable point exists 
to reach it. 
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